arXiv:1507.05943v2 [stat.AP] 26 May 2016 


MODELING THE PULSE SIGNAL BY WAVE-SHAPE FUNCTION 
AND ANALYZING BY SYNCHROSQUEEZING TRANSFORM 

HAU-TIENG WU\ HAN-KUEI WU^-^, CHUN-LI WANQ-^, YUEH-LUNG YANG®, 
WEN-HSIANG WU®, TUNG-HU TSAP AND HEN-HONG GHANG^’® 


Abstract. We apply the recently developed adaptive non-harmonic model based 
on the wave-shape function, as well as the time-frequency analysis tool called syn¬ 
chrosqueezing transform (SST) to model and analyze oscillatory physiological signals. 
To demonstrate how the model and algorithm work, we apply them to study the pulse 
wave signal. By extracting features called the spectral pulse signature, and based on 
functional regression, we characterize the hemodynamics from the radial pulse wave 
signals recorded by the sphygmomanometer. Analysis results suggest the potential 
of the proposed signal processing approach to extract health-related hemodynamics 
features. 


1. Introduction 

Hemodynamics is an essential ingredient of cardiovascular physiology which not only 
reflects the forces the heart needs to pump blood through the cardiovascular system, 
but also reflects the integrity of an individual’s physiological system. While there are 
many aspects to hemodynamics BElEiailS], it is common to evaluate hemodynamics 
by assessing the pulse wave signals recorded at different locations. For example, the 
brachial cuff systolic and diastolic blood pressures derived from the pulse wave have 
important clinical applications, long established beginning 1870 [6]. In addition to the 
blood pressure, the pulse wave signal itself contains an abundance of clinical informa¬ 
tion. For example, in subjects with dysfunctional left ventricle, different carotid pulse 
patterns can be observed, such as hyperkinetic pulse, pulsus alternans, pulse bisferiens, 
pulsus parvus et tardus, etc ma. □ However, it is difficult to obtain these carotid 
pulse wave signals non-invasively [U]. On the other hand, the non-invasive pulse diag¬ 
nosis (pulsology) based on, for example, the brachial or radial pulse signal, provides 
several aspects of physiological information, such as the central pressure wave mm- 
Other indices associated with the pulse wave signals recorded from different locations, 
including the augmentation index [5l Section 6.l|^ [H [12], are also available for clinical 
applications. Since there is a great amount of information within the pulse signal, a 
deeper and more extensive understanding of the pulse wave is undoubtedly important 
to better assess not only the cardiovascular but also the physiological integrity. 

^The hyperkinetic pulse is characterized by an increase in the velocity of the upstroke and amplitude. 
The pulsus alternans is a beat-to-beat variation in the amplitude of the pressure pulse. The pulsus 
bisferiens arterial pulse is perceived as two narrowly separated positive waves during systole. The 
pulsus pervious et tardus is a small and delayed arterial pulse. 

^The augmentation index indicates the incidence of reflected waves on the total pulse pressure. See 
[SJ Section 6.1] for a definition. 
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The common approaches to analyzing the pulse wave signals can be classified into two 
categories - time domain analysis and frequency domain analysis. In both categories, 
it is necessary to have either a representative pulse wave cycle, which reflects the 
interaction between the one time heart pump driving force and the impedance to 
blood flow, or a sequence of pulse wave cycles over a period of time, during which 
the heart rate is relatively constant. In the time domain analysis, different landmarks 
or features are identified from the pulse signal. For example, with the radial pulse 
signal is characterized by the height of the percussive wave, the height of the dicrotic 
notch, the length of the cardiac cycle, the length of acute ejection, etc. Different 
features with different physiological meanings can be obtained from the pulse signals 
recorded from other areas, and these features have been widely used in clinics for 
health evaluation [S]. Furthermore, new mathematical methods have been applied to 
investigate additional physiological information. For example, the central pulse wave 
could be reconstructed from the radial pulse wave by a generalized transfer function, 
and then be separated into inflection and reflection waves in order to evaluate markers 
of early vascular aging na. In the frequency domain analysis, the energy of different 
harmonic modes estimated by power spectrum analysis are also potential features for 
clinical applications. According to the resonance theory [H] , different harmonic modes 
are associated with the integrity of different organs [15]. Several existing works show 
the potential of this spectral approach, such as the relationship between the internal 
organ disorders and the spectral content [IS1[IZ|. Although this approach is commonly 
applied, it has two significant limitations. 

The first limitation is the over-simplification of the model underlying the analysis 
techniques. For example, heart rate variability (HRV) and pulsus alternans are of¬ 
ten not used in the analysis. Though the time domain landmark features reflect the 
hemodynamics from different aspects, they might be confounded by HRV and pul¬ 
sus alternans. Since HRV and pulsus alternans are inevitable physiological dynamics, 
without considering them in the model and analysis, both time domain and frequency 
domain analysis results may be inaccurate and the results might lead to incorrect inter¬ 
pretation. The second limitation is that it is not always straightforward to determine a 
representative pulse cycle or a period of time when the pulse wave signal is suitable for 
these analysis techniques, so that human intervention and subjective decision making 
are required. Indeed, to determine a pulse cycle, it is necessary to define what a pulse 
cycle means or, at least, to determine landmarks within the pattern to separate one 
pulse cycle from another. For example, with an electrocardiogram (ECG) signal, the 
repeating basic pattern is related to the electrophysiological activity of a normal heart 
beat, and a common chosen landmark for this is the R peak. As the R peak is usually 
dominant and significant, determining the oscillatory pattern for an ECG signal is not 
difficult. However, for pulse signal such as those acquired in this study, it is not always 
easy to define an oscillatory pattern or a landmark, in particular when the signal is 
recorded from an abnormal subject. Even more seriously, due to the possibility of a 
noisy record, the pulse cycle morphology can dramatically change from one pulse to 
another or from one subject to another (See Figure [^. Although this difficulty can be 
mitigated to some extent by noise reduction algorithms, due to the non-linear nature 
of the signal, their performance is not guaranteed. 
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Figure 1. An example of the radial pulse signal recorded from a subject 
with congestive heart failure. Note that at the 3’’'^ second there is an 
obvious artifact, indicated by the blue arrow, and at the 5.5*^ second 
there are two shock noise events. The relatively small pulse wave at the 
3'"'^ second is related to a premature atrial contraction. This is signihcant 
since patients with heart failure may have more premature beats than a 
normal control. 

To address these limitations, in this paper we apply a recently proposed and analyzed 
descriptive model, called the adaptive non-harmonic model, to describe the pulse wave 
signals. This model is characterized by what is referred to as the wave-shape function 
[18], which is dehned in order to capture what a pulse wave cycle should look like, 
including the features in the time domain and frequency domain mentioned above. 
The main characteristic of the proposed model is that the HRV and pulsus alternans 
are decoupled from the wave-shape function. A companion algorithm, referred to as the 
synchrosqueezing transform (SST), is applied to provide an accurate estimation of the 
wave-shape function. SST is a recently proposed time-frequency analysis technique that 
has a mathematical foundation and robustness to different kinds of noise and different 
practical issues have been well established in the literature uniEoj. By the adaptive 
non-harmonic model and SST, the above-mentioned limitations can be alleviated. By 
a suitable manipulation of the estimated pulse wave cycle, a vector can be obtained as 
an associated feature, which further leads to our final index, referred to as the global 
pulse signature (GPS) [T8l 1^ . 

This paper is organized as follows. In Section we systematically model the pulse 
signal using the notion of wave-shape function and the adaptive non-harmonic model. 
In Section we summarize the SST algorithm and discuss how the wave-shape func¬ 
tion can be estimated based on the SST and functional regression technique and then 
generate the GPS index. In Section we apply the proposed algorithm to study the 
radial pulse signal. The data collection procedure, the physiological background of 
pulsology and the analysis result are reported. In Section we discuss our hndings 
with their clinical interpretations and indicate directions for future research. 

2. Adaptive Non-Harmonic Model 

2.1. Adaptive non-harmonic model for the pulse signal. In this section, we 
propose a mathematical model to describe the pulse wave signal. It is well known 
that hemodynamics is a highly complex subject [21 El El El [23]. Different signals 
recorded from different areas of the body by different instruments can shed light on 
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different aspects of hemodynamics. The question of modeling hemodynamics has been 
discussed by several authors, who have proposed different models based on different 
physiological factors, such as pressure gradient, resonance, vibration, etc, have been 
proposed; see HI Ea El ID and the literatures cited therein for detailed discussions of 
these models. While these models shade a light on the understanding of the pulse 
wave signal, having a signal model to effectively quantify the pulse behavior remains 
a distant goal. In particular, inevitable physiological aspects, such as HRV and pulsus 
alternans, are commonly missed or disregarded in the existing models. Ignoring them 
in the model and analysis, however, might not only limit the analysis but also lead to 
incorrect hndings. In this paper, we apply a different model and approach to resolve 
these issues. 

The pulse wave signal is oscillatory in nature, and there are several well-known fea¬ 
tures that we consider here: how fast the signal oscillates, how much the oscillatory 
rate changes, how differently the signal behaves inside an oscillation, etc. Understand¬ 
ing these features plays a signihcant role in the data analysis, which can help to obtain 
more information about the physiological system. Based on these observations, we con¬ 
sider the following adaptive non-harmonic model [18] to model the pulse signal. This 
model differs from most existing models in that it is purely phenomenological] that 
is, the parameters in the model are solely driven by observations of the physiological 
signal, rather than by explicit and quantitative underlying mechanisms. 

2.1.1. Wave shape function. Based on the periodic nature of an oscillatory function, 
we introduce the wave shape function, which is used to model how a signal oscillates 
over the period. This idea has been previously studied in [THl ES] ES] . To resolve the 
difficulty of defining a period, we shall call a function / t — periodic, where r > 0, if (1) 
the periodic condition is satisfied; that is, for all t G M and all fc G Z, f{t + kr) = f{t)] 
and (2) for all 0 < r' < r the periodic condition is not satisfied. Take two positive 
values 6 and 6 and a positive integer D. A 1-periodic function s is called a wave shape 
function with parameters S, D and 6 if the following conditions are satisfied. First, 
the function is differentiable and its derivative is a Holder continuous function with 
the Holder coefficient a > 1/2. Denote 's to be the Fourier transform of the function 
s. The second condition that a wave shape function should satisfy is that s has the 
unit energy (unit L^-norm) and all the Fourier modes 's{k), k ^ 1 are dominated by 
the product of 6 and the hrst mode coefficient; that is, 

(1) V/c G N, with k ^ 1, |s(/c)| < 5 |s(l)| . 

Furthermore, it is desirable for s' to be mostly concentrated in the low frequency region, 
which is quantihed by 9 by: 

(2) ^ |n?(n)| < 9. 

n>D 

Conditions ([^ and (|^ reflect the practical finding of the spectral analysis of the pulse 
signals [23 EH El] , and it can be observed that the amplitudes of tenth and higher or¬ 
der harmonics are negligible. Note that the commonly selected landmarks and lengths 
considered in the time domain analysis could be understood as morphological features 
describing the wave-shape function that is modeling how a pulse repeats itself. While 
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there is a one-to-one relationship between these landmarks and Fonrier coefficients 
's, this relationship is nonlinear in nature. Thus, while we could view the Fourier 
coefficients of the wave shape function as the “features” of the pulse wave, their physi¬ 
ological interpretations are not directly related to the hemodynamic interpretations of 
these landmarks. See Figure]^ for an example of the wave-shape function describing 
the pulse wave signal. 

2.1.2. Adaptive non-harmonic model. Using the notion of wave shape function, we 
now describe our phenomenological model to capture the recorded pulse signal. Fix 
parameters <5, 9, D for a wave shape function s. We consider the following Intrinsic 
Mode Type (IMT) functions to model the pulse wave signal. An IMF function, /, 
is a bounded and continuous function with a continuous derivative that satishes the 
following format: 

(3) f{t) = A{t)s{2TT(j){t)), 

where A is a positive differentiable function and 0 is a monotonically differentiable 
function. Intuitively, A{t) describes how large the oscillation is at time t, and the 
positive function 4>'{t), the hrst derivative of 0, describes how fast the oscillation is at 
time t. To see how it is interpreted in this way, consider a constant positive function 
A and a linear function with a positive slope as 0. Suppose 0(f) = f^t, where > 0- 
In this case, we know that / is a harmonic function with the frequency cind the 
amplitude A. We could thus view (|^ as a generalization of the harmonic function. 
Though the heart rate is not constant, normally it does not change suddenly. Therefore 
we need following conditions to better quantify the pulse signal. Fix a small positive 
constant e. Then we assume that A{t) does not change too fast in the sense that its 
derivative is bounded by e0'(f) and that 0'(f) does not change too fast in the sense 
that 0"(f) exists and is bounded by e0'(f); that is, we have 

(4) ^ and |0^^(f)| < 60'(f) for all time t. 

We would thus model the pulse wave signal by ([^ with the condition (|^, and call this 
model the adaptive non-harmonic model. 

We would call the monotonically increasing function 0(f) the phase function, fhe 
instantaneous frequency (IF), and A(f) the amplitude modulation (AM). An important 
issue regarding the identihability issue of the phase function, the IF and the AM 
cannot be ignored if the discussion is to be rigorous. Indeed, there might be inhnitely 
many different ways to represent a cosine function go{t) = cos(27rf) in the format 
a(f) cos(27r0(f)) so that a > 0 and 0' > 0, even though it is well known that g^it) is a 
harmonic function with amplitude 1 and frequency 1. Indeed, we could hnd inhnitely 
many smooth functions a and 0 so that g^it) = cos(f) = (1 -f- Q;(f)) cos(27r(f -|- 0(f))), 
and in general there is no reason to favor a{t) = /3{t) = 0. Clearly, before resolving 
this issue, amplitude 1 and frequency 1 cannot be taken as reliable features to quantify 
the signal go. This identihability issue has been well studied in [201 EQ], hnding that if 
go{t) = A(f)s(0(f)) and go{t) = (A(f) -f a;(f))s(0(f) -|- 0(f)) also holds, and both satisfy 
the condition (|^, then |a(f)| and |0'(f)| are both bounded by e for all time f G M. 
Note that the IF and AM are always positive, but usually not constant. Clearly, 
when 0 is a linear function with a positive slope and a is a positive constant, then 
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the model is reduced to the harmonic model and the IF is equivalent to the notion 
frequency in the ordinary Fourier transform sense. The conditions |aXt)| < and 

W'it)\ < force the signal to locally oscillate “regularly”, that is, around time to, 

the signal A{t)s{((){t)) oscillates like A{to)s{((){to) — + 4>'{'to)t), and hence the 

nominations of IF and AM. We mention that this model is a special case of a wider 
class of model composed of multiple oscillatory components in the sense that we only 
have one oscillatory component in the model. For details about a more general model, 
see [mUHlEOlED]. 

2.1.3. Physiological interpretation. The main reason we consider a time-varying fre¬ 
quency and amplitude model like ([^ is to capture the physiological facts of HRV and 
pulsus alternans. 

Since the heart does not beat at a constant rate [21], the pulse signal should not 
oscillate with a constant frequency. The non-constant heart rate is modeled by IF 
and hence the HRV could be evaluated from analyzing IF. In the past decade, due to 
the health information embedding in HRV and the trend towards personalized health 
care requirements, estimating the HRV from the pulse signal extracted from different 
resources has attracted signihcant research [32[ |33l |32] . However, the existence of HRV 
is commonly ignored in the pulse analysis literature. For example, in the traditional 
spectral analysis approach to the radial pulse, researchers need to analyze the pulse 
signal over the interval where the heart rate closely resembles a constant [2HI US] . 

The possible pulsus alternans is captured by the AM. This model is particularly 
important in patients susceptible to heart failure since pulsus alternans is a manifes¬ 
tation of severe impairment of the left ventricular systolic function [7] . An example of 
a wave-shape function and an IMT function representing the pulse signal are shown in 
Figure 



Figure 2. Left: an illustrative wave-shape function s with D = 5, 

6 = 0.59 and 6 = 0. Right: an IMT function with the wave-shape 
function s and time varying amplitude and frequency. Here the frequency 
is about 1.2 Hz. 

2.1.4. Recorded pulse signal. In practice, noise is inevitable when recording signals. 
Thus, we model the recorded pulse signal by 

(5) V(t) = A(t)s(0(t))+ cr(t)$(t), 

where A{t)s{(j){t)) is the adaptive non-hamornic model for the pulse signal, cr is a 
slowly varying smooth function quantifying the possible non-stationarity in the data 
collection process, and $ is a stationary random process with unit standard deviation 
describing the noise. See, for example, [20] for a discussion of the noise. Here, we do not 
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assume that the noise is Gaussian or white, while the commonly encountered Gaussian 
white noise is a specihc example. Since in general the noise might be complicated with 
time-dependence and be far from Gaussian, we take this noise model into account. 

While the pulse signal could be recorded by different ways, like tonometer, photo- 
plethymography, video [53], etc, the wave-shape function could be different for different 
types of pulse signals. 


3. Methodology 

3.1. Synchrosqueezing transform. It is a common practice to apply the Fourier 
transform to study oscillatory signals, in particular the pulse wave signal in which 
we are interested unmsi Ej. As useful as the spectral analysis is, however, it is 
well known that when the signal is not composed of harmonic functions, the power 
spectrum determined by the Fourier transform might be misleading. For spectral pulse 
wave analysis, since HRV and pulsus alternans are inevitable, in practice analysts must 
carefully choose the signal. But there is still no guarantee that the HRV and pulsus 
alternans influence could be eliminated, and this could become a confounder in the 
analysis. 

While the oscillatory signals with time-varying AM and IF are ubiquitous, much 
effort has been extended over the past few decades to address this problem, and several 
approaches have been proposed, ranging from empirical mode decomposition and its 
variations [351 ESI ED, to optimized window and the approximation theory approach 
[38l [39] . to name but a few. Moreover, research interest has recently been extended 
to multivariate time series analysis to further take the spatial 

information into account. Among the different approaches, one active field regarding 
this issue is time-frequency (TF) analysis. Based on different principles, several TF 
analysis techniques have been proposed. Typical examples include linear methods such 
as the short time Fourier transform (STFT) and the continuous wavelet transform 
(GWT), or the quadratic methods such as the Wigner-Ville distribution or the Gohen 
class. We refer the reader to a standard textbook m for more information. While 
these methods have attracted a great deal of attention in different helds other than 
signal processing, they are limited by either the Heisenberg uncertainty principal or the 
mixing issue. We demonstrate this idea by discussing STFT. The main idea forming the 
basis for STFT is dividing the signal into overlapping small pieces, and then studying 
the spectral behaviors over these small pieces. Mathematically, for a function /, its 
STFT associated with a window function h{t) can be dehned as 

Vj^\t,r]) ■= j /(s)/i(t - s)e-*2-^(*-")ds 

where f G M is the time, r] G M"*" is the frequency, h is the window function determined 
by the user - for which a common choice is the Gaussian function with kernel bandwidth 
cr > 0, i.e. h{t) = (27rcT)“^/^e“* . However, the Heisenberg uncertainty principal 

limits how well the spectrum could be estimated over each piece, thereby limiting the 
STFT. Similar discussions hold for GWT and other linear TF analysis techniques, and 
we refer the reader to [13] for details. In the current work, in order to capture the 
hemodynamics, which have oscillations on the order of 1 second, we run STFT with 
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a = 0.5 so that the window is not too long to lose the dynamic information, and is not 
too short to cause a numerical artifact. 

To handle these fundamental issues, a nonlinear TF analysis technique, synchrosqueez¬ 
ing transform (SST) [T^l ED], which is a special reassignment technique (RM) |1D], is 
proposed to obtain a sharper time-frequency representation. Since these techniques 
were introduced, they have been widely applied in different helds; see 071 for a sum¬ 
mary of its applications in different fields. 

Here we briefly summarize SST, and refer the reader to [16| for RM. We mention that 
the SST could be applied to different linear TF analysis, like STFT, CWT, wave packet 
transform [25] or S-transform ^8] and theoretically the results do not depend signifi¬ 
cantly on the chosen linear TF analysis. But to simplify the discussion, we choose to do 
the analysis with STFT. The Matlab code for the SST algorithm based on both CWT 
and STFT can be downloaded from 
home/download. In brief, the SST 
STFT by reallocating the STFT coefficients along the frequency axis to the “correct” 
frequency slot which represents the IF of an oscillatory component. Mathematically, 
The SST of / is dehned as 


https://sites.google.com/site/hautiengwu/ 
sharpens the TF representation determined by 


Sf\t,0 ■= lim 

^ rv — 


9a{\^ - dr], 


where ga is an approximate ^-function in the sense that g is a fast decaying smooth 
function with / 5 f(a:)da: = 1, so that ga(t) := ^g{^) tends to the Dirac delta measure 6 
supported at 0 weakly as a —)■ 0, and with dehned by 



— OO 


if ^ 0 

otherwise. 


In plain language, by reading we collect all STFT coefficients indicating 

the existence of an oscillatory component with frequency ^ to the slot Note that 
compared with RM, in SST, the coefficients are reallocated along the frequency axis, 
so the causality is preserved; second, in SST we reallocate the STFT coefficient instead 
of the spectrogram coefficient. These two facts allow us to reconstruct the oscillatory 
components of interest, in particular when the signal is noisy or is composed of several 
oscillatory components. It has been clearly found that for the phenomenological model 
of the pulse wave signal of interest, the IF 4>\t) and the AM A{t) can be accurately 
estimated from the recorded pulse signal HSlEo]. Precisely, we could prove that at 
time t, the coefficients in are dominant when ^ This property allows 

us an accurate estimate of the IF 0' by, for example, the curve extraction technique. 
Denote the estimated 0' by 0'. We can then estimate the amplitude modulation A{t) 
and the phase function 0(f) by building 

(6) R{t) -h{0)-' [ _ 

The estimator of A{t) is thus dehned as A{t) := |R(f)|, and hence an estimator for 0(f), 
denoted as 0(f), can be obtained by unwrapping the phase of the complex-valued signal 





WAVE-SHAPE FUNCTIONS AND PULSE SIGNALS 


9 


R(t)/A(t). We refer readers interested in SST to [191 120] for the detailed numerical 
algorithms and the theory beyond them. 

Here, we show the results of pulse analysis by applying the SST in Figure To 
demonstrate its beneht on a noisy signal, we artihcially add noise on the second half of 
the signal, and show the result in Figure]^ In this example, the noise is an ARMA(1,1), 
where ARMA stands for autoregressive and moving averaging, time series determined 
by the autoregression polynomial a{z) = 0.5z +1 and the moving averaging polynomial 
b{z) = —0.3z+l, with the innovation process taken as i.i.d. student random variables 
so that the signal to noise ratio, dehned as 20 fog is 0 dB. It is clear that though 

the signal to noise is low, the noise is non-st at ionary since it exists only over a hnite 
interval, and the noise has the fat tail behavior described by the student random 
variables. Thus, the SST algorithm could reliably extract the instantaneous frequency, 
so that we could obtain reliable HRV information. 


3.2. Feature extraction by estimation of the wave shape functions. As dis¬ 
cussed in Section |2.1[ the main feature of interest for the pulse analysis is the wave¬ 


shape function in the adaptive non-harmonic model, in particular in the pulse spectral 
analysis. In this section, we describe an algorithm to extract the wave shape function. 
Note that since the wave-shape function can be expanded by its Fourier coefficients, 
the pulse signal f{t) = A{t)s{cj){t)) can be represented as 


(7) f{t) = A{t) ^ [ai cos(27r70(f)) + /3i sin(27r£0(f))], 

where G M and G M are the Fourier coefficients of the shape function s. In this 
study, we estimate the wave shape functions based on Q using the standard functional 
regression HHIEI]. 

Consider the wave shape functions s with parameters 5, id, 6^ in ([^. To simplify our 
discussion, we assume that 6 = 0 and that the noise is stationary, that is a = 1. We 
would choose Zi = 6 in this study, based on the discussion in Section Thus, the 
pulse signal Q becomes 

D 

f{t) = A{t)ao -|- A{t) [«£ cos(27ri'0(t)) -f sin(27r70(f))]. 

1=1 

After discretization by the sampling period At > 0 over time interval [At, NAt], the 
recorded pulse signal is saved digitally as a A-dim vector, Y G so that its /-th 
entry is 1^/ = /(/At) -|- ^i, where / = 1,... ,N and $ is a random vector satisfying 
var(<I);) = 1 for all /, which might not be Gaussian and the covariant matrix might not 
be the identity. Denote the discretized estimators of A{t) and /(t) by A G and 
(j) G Note that it has been well established the discretized estimation of A and 0 
are accurate with error of order e m- To simplify the discussion, we assume that the 
estimates A and (j) are precise without error; that is, A{1) = A{lAt) and /(/) = /(/At) 
for all / = 1,..., A. In the general case, the analysis result will deviate by an error of 
order of e. We thus construct the following “functional vectors” 


c=[clcl...,cl,dl...,direR^^-^ 
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Figure 3. (a) The radial pulse signal R{t) recorded from a normal 
subject. It is also clear that the IF and amplitude modulation are not 
constant due to the inevitable HRV and pulsus alternans. The peaks 
are marked by gray circles, and the lengths between two consecutive 
peaks are shown. Note that there are two peaks around the second 
second, and we choose the one with the larger value as the peak, (b) 
The time-frequency (TF) representation of R{t) determined by the syn¬ 
chrosqueezing transform. The dominant curve in the TF representation 
is associated with the IF induced by the heart rate variability (HRV). 
It is clear that the pulse around the 5*^ second takes a longer time to 
hnish, which leads to the slower instantaneous frequency. Indeed, the 
y-axis of the dominant curve around time 5.3-th second is 1/1.15 = 0.87, 
which reflects how fast the signal oscillates at that moment. Note that 
the artifacts around the 5.5 second and the 7.5 second do not play a 
major role in the analysis results. 


where £ = 0,..., H and q, are iV-dim vector whose k-ih. entries are 
Ciik) := A{k) cos(27r£0(A;)), d^ik) := A(k) sin(27r£0(A;)), 
where /c = 1 ,..., iV. As a result, the recorded pulse signal satishes 

(8) V = 7^c + $ 

where 7 = [dq, Di, ..., ctd, A, • • •, fdo] e 

To estimate the parameters and /3i from the functional vectors c, observe that 
the {2D -|- 1) X {2D + 1) matrix cc^ is diagonal dominant. Since E(At<l>c^) = 0 and 
var(At<hc'^) = 0{LAt). Thus we can estimate 7 by 

(9) j := {Yc^){cc^)-\ 
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Time (Second) 

Figure 4. (a) The radial pulse signal R{t) recorded from a normal 
subject that is contaminated by the autoregressive and moving averaged 
noise generated from the student ts i.i.d. random process from the 2 .5-th 
second to the 5.5-th second, which is denoted as Y{t). The radial pulse 
signal is the same as that shown in Figure Note that there are several 
spikes in Y[t), which are generated by the fat tail natural of the student 
t 3 random variable, (b) The time-frequency (TF) representation of Y[t) 
is determined by the synchrosqueezing transform (SST). The dominant 
curve in the TF representation is associated with the IF induced by the 
heart rate variability (HRV). It is clear that the TF representation deter¬ 
mined by SST is robust to the noise even if the noise is time-dependent 
with a fat tail distribution, while the dominant curve representing the 
instantaneous frequency is slightly deformed due to the noise. 

where 7 = [5o, 5i,..., 5 d, Pi, ■ ■ ■, Pd]"^ £ and hence estimate the oscillatory 

signal by 

f{t) ■= 7^c. 

When / contains an accurate estimation of the wave-shape function s with error of 
order e, we could take the Fourier coefficients 7 into account as the feature of the pulse 
signal. We call the 2D + 1 dimensional vector 7 the spectral pulse signature (SPS) 
for the recorded pulse signal. Note that the £-th component of the power spectrum 
of s could be estimated by (5^ - 1 - P‘f)/Y The main benefit of this approach is that 
the influence of HRV, which is modeled by IF, is eliminated automatically and the 
wave-shape function is better estimated. 

3.3. Comparing the present method to previous ones. Here we summarize the 
difference between our approach and the commonly applied spectral analysis mm 
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[29]. Recall that in spectral analysis, the power spectrum of a selected recorded pulse 
signal is evaluated and the energy of different harmonic modes are considered as fea¬ 
tures of the subject, indicating aspects of physiological health. The selection criteria 
for the pulse signal interval is that the heart rate is almost a constant [291IT9] , since the 
power spectrum approach is sensitive to the non-constant frequency and non-constant 
amplitude signal. The hrst difference is that based on the adaptive non-harmonic 
model and SST, the instantaneous frequency and amplitude modulation can be accu¬ 
rately obtained, so we do not need to select an interval from a recorded pulse signal. 
Second, since the HRV and pulsus alternans are physiologically inevitable, the proposed 
approach is more physiologically feasible. 

Third, note that the information obtained from the spectral analysis is different from 
the SPS. Indeed, under the assumption for (|^, we have the following direct expansion: 

D D 

s{t) = ^ = ao + ^ [«£ cos(27r£f) -1- f3e sm{27iit )], 

£=-D e=i 

where ao = s(0) and for £ = 1, ..., D, s{£) = ae > 0, G [0, 27r), ae = 2a£ cos(6e) 

and /3e = —2aesm6£. Here, 6^ reflects the phase of the £-th harmonic component 
hidden inside the wave shape signal. Clearly, since power spectrum analysis commonly 
takes |s(f')p = aj, the phase information of the Cth harmonic component is missed. 
However, in SPS the phase information is preserved and the hemodynamics could be 
more faithfully captured. 

4. Testbed: the radial pulse diagnosis 

To test how well the proposed adaptive non-harmonic model and the algorithm work 
in practice, in this section we study the radial pulse wave signal on congestive heart 
failure subjects. The signals are recorded from three landmarks. First, the radial 
styloid process; second the middle position between the radial styloid process and 
the palm and the proximal point with the same distance from the hrst to the second 
location. We mention that these three locations are commonly recognized in the pulse 
diagnosis; the hrst one is called guan, the second one is called chun and the third one 
is called chi. We thus use guan, chun and chi to refer to these three locations in this 
study. 

4.1. Material. All protocols in this study were approved by the Institutional Review 
Board of Chang Gung Memorial Hospital, Linkou, Taoyuan, Taiwan (93-6288), Taiwan 
and written informed consent was obtained from all patients. Nineteen normal subjects 
without history of heart disease are included in the control group, and 17 subjects 
with congestive heart failure (CHF) are included in the study group. The diagnosis of 
CHF subjects is based on the criteria indicated by the Framingham heart study. The 
participants were invited for pulse examination in a quiet and temperature-controlled 
room in Chang Gung Memorial Hospital, Linkou branch, Taiwan. Pulse wave signals 
were recorded from chun, guan and chi positions of both hands by a tonometer (Wang’s 
sphygmometer, PDS-2000). The sampling rate of the signal is at lOOHz. For each 

^In the literature sometimes chun is called inch or chon, guan is called bar or gwan and chi is called 
cubit or cheok. 
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subject, we collect 10 seconds signal for each position on both hands, and repeat 2 or 
3 times. The pulse wave was recorded in sitting position with the wrist comfortably 
resting on a small pillow at the level of the heart 

We recruited 17 patients with CHF for the study group and 19 normal individuals 
for the controls. The age of the study and control group are 64.3 ± 23.7 and 63.2 ± 15.8 
respectively. The male/female ratio were 15/7 in the study group and 10/10 in the 
control group. There were no signihcant differences in age or sex. As a result, we 
obtain 53 (respectively 39) pulse wave signals recorded from chun from the normal 
(respectively CHF) group, 55 (respectively 34) pulse wave signals recorded from guan 
from the normal (respectively CHF) group and 50 (respectively 41) pulse wave signals 
recorded from chi from the normal (respectively CHF) group. 

4.2. Statistical Analysis and Global pulse signature (GPS). To test if the SPS 

indices of the normal group differ from those of the CHF group, we apply the cur¬ 
rently developed one-way ANOVA for functional data, called the globalized pointwise 
F (GPF) test [50]. For readers having interest in this technique, we refer them to [50] . 
In this study, we consider p values less than 0.01 as signihcant. 

As the SPS index is a high dimensional vector, it is not easy to visualize. To provide 
an easy-to-use index for the pulse diagnosis, we consider the following approach to 
integrate the information in the SPS index. Suppose from a hxed position of the i- 
th subject we obtain a SPS index 7 ^. The associated outcome of this SPS index is 
denoted as yi, and yi = I (respectively yt = 0) means the subject is in the CHF group 
(respectively control group). Thus we have the dataset {ji,yi}^i. When the SPS 
index is located in the 2D -I- 1 Euclidean space and the sampling size is limited, we 
apply the partial least squares (PLS) regression to hnd a linear regression model by 
projecting SPS and the response variables to a new space. Here we briehy recall the PLS 
regression. PLS regression Ends components from X = that are relevant for 

the outcome 3 ^ = {yi}f^i by seeking a set of components that performs a simultaneous 
decomposition X and y with the constraint that these components explain as much as 
possible of the covariance between X and 3^. Then the decomposition of X is applied 
to predict the group. For details about PLS regression, see [51]. Suppose the PLS 
regression coefficient is /3 G Then the prediction result under the linear 

regression model for 7 j, denoted as iji := [1 7 J /3 G M, is referred to as the global pulse 
signature (GPS) index. The GPS index is integrated information derived from the SPS 
index, which rehects the subject’s condition of interest. 

For the purpose of prediction based on GPS, we can further apply the receiver 
operating characteristic (ROG) to determine the threshold to classify the subjects into 
two groups. We report the sensitivity, specihcity, accuracy and the area under curve 
(AUG) to evaluate the classihcation result. The conhdence interval (GI) of AUG is 
evaluated by 1000 bootstrap replicas. To assess how the results based on PLS and 
ROG will generalize to an independent data set, we run leave-one-out cross validation 
(LOOGV) 200 times for each position of both hands, and report the accuracy. 

4.3. Results. First, we show the synchrosqueezing transform of the pulse wave signal 
from a subject with GHF and the associated estimated wave shape function in Figure 

Note that due to the inevitable deviation, for example the one at the 4-th second. 
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and the HRV, the power spectrum estimated from the pulse wave signal is spreading. 
Note that estimating the wave shape function could be viewed as the power spectrum 
analysis of the pulse wave signal after correcting the IF and AM. In other words, as we 
could estimate IF and AM accurately from the pulse wave signal, we could resample 
the signal according to the estimated IF and then normalize the signal by the estimated 
AM. 


(a) 



(b) 




(d) 


, . 


1 2 3 4 5 6 

Frequency (Hz) 


Figure 5. (a) The pulse signal recorded from a subject with conges¬ 
tive heart failure (CHF), with the y-axis as the arbitrary unit, (b) 

The estimated wave shape function for the subject with CHF. (c) The 
time-frequency (TF) representation determined by the synchrosqueezing 
transform. It is clear that the pulse around the 2.5*^ second takes a longer 
time to hnish, which leads to the slower instantaneous frequency (the 
dominant curve in the TF representation is around 0.9 Hz at time2.5*^ 
second). It is also clear that the instantaneous frequency is not con¬ 
stant due to the inevitable heart rate variability (HRV). Note that while 
there is a signihcant deviation at the 2.5-th second, the estimation result 
catches most of the shape information, (d) The power spectrum of the 
estimated wave shape function shown in (b). 

The set of SPS indices of different groups from different positions is shown in Figurej^ 
We could see that the means of the normal group and the CHF group are different. The 
GPF functional ANOVA test shows that the SPS indices evaluated from all positions 
on both hands are signihcantly different. Except for the p value of the chun position 
on the right hand, which is 4.4 x 10“'^, the p values of other positions are < 10“^. 

Then we apply PLS to obtain the GPS index to distinguish the two groups. The 
histogram of GPS indices determined from different positions on both hands are shown 
in Figure It is clear that the GPS of subjects in the GHF group is smaller than that 
of the control group. The ROG analysis results, including the sensitivity and specihcity 
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Figure 6 . The deviation from the mean of each entry of the spectral 
pulse signature (SPS) is shown in the error bar plot. The error bar is 
determined by the standard deviation. In the top row, the distribution 
of SPS determined from guan, chun and chi positions of the left hand 
are shown from left to right. In the bottom row, the distribution of SPS 
determined from guan, chun and chi of the right hand are shown from left 
to right. The results of the normal subjects are shown in the gray error 
bar, and those of the subjects with congestive heart failure are shown in 
the black error bar. 


and AUC, from different positions are shown in Figure]^ The sensitivity, specificity, 
accuracy, AUC and the accuracy of LOOCV of the GPS determined from different 
positions on both hands are summarized in Table 



GPS GPS GPS 

Figure 7. The histograms of GPSs determined from chun, guan and 
chi positions on the left (respectively right) hand are shown on the left, 
middle and right subfigures on the top (respectively bottom) row. The 
normal group is shown in the gray color and the congestive heart failure 
group is shown in black. Gray stars represent the determined thresholds 
by the ROC binary classihcation. 
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Figure 8 . Left (respectively middle and right): the receiver operating 
cnrve (ROC) of the global pnlse signatnre (GPS) of the chnn (respectively 
gnan and chi) position on the left hand is shown in the black cnrve and 
that of the right hand is shown in the gray hand. Black and gray circles 
are the optimal operating point of the ROC. 



Left hand 

SEN 

SPE 

AC 

AUC 

LOOCV 

chnn 

0.87 

0.9 

0.88 

0.92 (0.81-0.96) 

0.8 

gnan 

0.89 

0.79 

0.85 

0.91 (0.85-0.96) 

0.74 

chi 

0.94 

0.78 

0.87 

0.89 (0.8-0.95) 

0.74 


Right hand 

SEN 

SPE 

AC 

AUC 

LOOCV 

chnn 

0.87 

0.79 

0.84 

0.89 (0.79-0.95) 

0.74 

gnan 

0.94 

0.73 

0.85 

0.91 (0.83-0.95) 

0.79 

chi 

0.94 

1 

0.96 

1 (0.97-1) 

0.94 

. Sensitivity 

(SEN), speci 

icity (SPE), accuracy (AC), 


Table 1. 

nnder cnrve (AUC) and the accnracy of the leave-one-ont cross validation 
(LOOCV) of the global pnlse signatnre evalnated from different positions 
on both hands. The confidence intervals of the AUC are reported in 
parentheses. 


5. Discussion 

In snmmary, in this paper we stndy the physiological signal by the adaptive non¬ 
harmonic model, the SST and the fnnctional regression techniqne. The nsefnlness 
of the proposed scheme is snpported by an enconraging analysis resnlt of the radial 
pnlse signal. Althongh we analyze the radial pnlse signal as the test case for this 
stndy, it shonld be noted that the proposed model and analysis techniqne could be 
applied to other pulse signals obtained by different instruments. For example, con¬ 
tact photoplethysmogram (PPG) measurement [52] or PhysioCam non-contact PPG 
measurement [33] , which represent the changes of blood volume in the vessel obtained 
through an optical transmission measurement or real-time camera images, and hence 
reflect the change in vascular blood volume associated with the cardiac beat. While 
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these signals represent different aspects of the hemodynamics than the radial pulse 
signal we study in this paper, we could expect to obtain a broader angle of view about 
human health if information obtained from these signals could be combined. We will 
report the research progress in the future work. Since the potential of the SST and 
other nonlinear TF analysis techniques have been shown in this study and other clinical 
problems, for example, ISSlEllES], by taking the wave-shape model into account, we 
could expect a broader application and better analysis results. 

We also point out the relationship of the pulse signal analysis results with the tradi¬ 
tional Chinese medicine (TCM) theory. In short, the results we obtained in this study 
could lead to a potential means to help establish the foundation of TCM theory. See 
Figure for an illustration of the summary of the pulse diagnosis theory in the TCM 
theory. Note that while TCM has been widely applied in the eastern culture, up to now. 



Figure 9. An illustration of anatomical locations associated with the 
terminologies used in pulse diagnosis. In the TCM pulse diagnosis, it 
is stated that the pulse on the right hand manifests the condition of Qi 
[SniETj, while the pulse on the left hand reflects the blood. For different 
pulse positions, left chun, guan and chi reflect the heart, the liver, and 
the kidney respectively; the right chun, guan and chi reflect the lung, the 
spleen, the kidney respectively. The left chi manifests mainly the kidney 
yin, and the right chi manifests mainly the kidney yang (the Life Gate 
in TCM). 

a systematical/scientific theory understanding the mechanism beyond pulse diagnosis 
is not yet well established [SB], but its usefulness has been demonstrated [SB]. There 
are several studies based on hemodynamics aiming to understand the mechanism; for 
example, Wang et. al. PI [59] proposed that the pulse wave consists of numerous 
harmonic waves, and each harmonic wave is associated with an internal organ and car¬ 
ries information of different meridians over the body. That study also noted that the 
harmonic waves of the upper, middle and lower section of the body may correspond to 
the three sections of pulse in timeline and that the waveform on both hands are not 
totally the same. On the other hand, our approach is purely from the phenomenologi¬ 
cal viewpoint and adaptive harmonic analysis. In our results, the classification result 
of the CHF group and the normal group is better on the right hand side, especially 
the right chi. Since the right chi position demonstrates the most significant difference 
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between groups, it may reveal the importance of the right chi (kidney yang) signal 
in evaluating the clinical condition of CHF. According to TCM theory, kidney yang 
is the fundamental support to maintain the human life. In Western medicine, the 
hemodynamics of renal artery, such as the blood pressure, plays an important role in 
regulating the overall whole cardiovascular function. On the other hand, a possible 
cause for this observation is that the pulse pressure of the right arm is usually higher 
than the pulse pressure of the left arm. It is presumed that the pulse wave signal of 
the right arm may have a higher signal-to-noise ratio than the pulse wave signal of the 
left arm. The left chun and right guan position are the two second important roles in 
our analysis. According to TCM theory, the left chun position (heart) manifests the 
general condition of the cardiovascular system, hence our results are consistent with 
clinical experience. Thus, our hndings partially support the pulse diagnosis theory in 
TCM that the waveforms on different positions of radial artery contain different infor¬ 
mation. A similar hnding was reported by Young et. ah, who found that the although 
the augmentation indices determined from the radial pulse waves recorded from chun, 
guan and chi are not signihcantly different, the estimated aortic augmentation indices 
determined from the radial pulse waves recorded from chun, guan and chi are not iden¬ 
tical [60]. In conclusion, since the study of pulse diagnosis from different aspects is an 
active held, we would expect our proposed model and method could help to further 
study the experiences and working practice of pulse diagnosis, e.g. the nature and 
dynamic of disease, and its relationship to modern hemodynamics. 

Limitations of this study should also be mentioned. First, the tonometer (PDS-2000) 
we applied in this study records only the two-dimensional data (pressure-time) of the 
pulse. Since more advanced instruments now available can obtain three dimensional 
data [61], further study should be undertaken. Second, although the phenomenological 
model we propose is capable of capturing IF and AM as well as the wave shape function, 
it is clearly not the optimal solution. It is clear that there is a more complex interaction 
between IF, AM and the wave shape function than what we consider in the adaptive 
non-harmonic model. On the one hand, a more general model, like the time-varying 
wave shape function, could be considered based on the physiological background. On 
the other hand, we conjecture that this relationship might be better captured by com¬ 
bining the existing hemodynamic models with the proposed phenomenological model. 
This hner model might better capture the physiological information hidden inside the 
pulse wave signal and lead to a better algorithm to better study the recorded pulse 
wave signal. A systematic study and its application of this issue will be reported in 
future work. Moreover, from the clinical viewpoint, the sample size in this study is 
limited, and the interesting clinical problems, like outcome prediction or early CHF 
diagnosis are not discussed. Also, to simplify the discussion and avoid possible con- 
founders, in this study we limit our analysis to subjects with CHF. Thus, we could 
not make the hnal conclusion about the pulse diagnosis. To use the research results in 
clinics, a larger scale clinical study with CHF and other diseases is needed to conclude 
the current hndings, and we will report our continuing research in the future. 
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